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Abstract 


We present constraints on cosmological parameters from the Pantheon+ analysis of 1701 light curves of 1550 
distinct Type Ia supernovae (SNe Ia) ranging in redshift from z = 0.001 to 2.26. This work features an increased 
sample size from the addition of multiple cross-calibrated photometric systems of SNe covering an increased 
redshift span, and improved treatments of systematic uncertainties in comparison to the original Pantheon 
analysis, which together result in a factor of 2 improvement in cosmological constraining power. For a flat 
ACDM model, we find Qy=0.334+40.018 from SNe Ia alone. For a flat woCDM model, we measure 
Wo = —0.90+0.14 from SNela alone, Hyp =73.5+1.1kms 'Mpc~' when including the Cepheid host 
distances and covariance (SHOES), and wo =—0.978*903} when combining the SN likelihood with Planck 
constraints from the cosmic microwave background (CMB) and baryon acoustic oscillations (BAO); both wo 
values are consistent with a cosmological constant. We also present the most precise measurements to date on 
the evolution of dark energy in a flat wow,CDM universe, and measure w, = =O. from Pantheon+ SNe Ila 
alone, Hy = 73.3+1.1kms 'Mpc ! when including SHOES Cepheid distances, and w,=—0.65*)35 when 
combining Pantheon+ SNe Ia with CMB and BAO data. Finally, we find that systematic uncertainties in the use 
of SNe Ia along the distance ladder comprise less than one-third of the total uncertainty in the measurement of Hp 
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and cannot explain the present “Hubble tension” between local measurements and early universe predictions 


from the cosmological model. 


Unified Astronomy Thesaurus concepts: Cosmology (343); Dark energy (351); Dark matter (353); Type Ia 
supernovae (1728); Cosmological models (337); Expanding universe (502) 


Supporting material: machine-readable table 


1. Introduction 


Type Ia supernovae (SNe Ia) anchor the standard model of 
cosmology with their unmatched ability to map the past 10 
billion years of expansion history. SNe Ia provided the first 
evidence of the accelerating expansion of the universe (Riess 
et al. 1998; Perlmutter et al. 1999), and they remain invaluable 
because they are (1) bright enough to be seen at large cosmic 
distances, (2) common enough to be found in large numbers, 
and (3) can be standardized to ~0.1 mag precision in brightness 
or ~5% in distance per object. 

Statistical leverage from large samples of SNe Ia has grown 
rapidly over the last three decades, and well-calibrated and 
standardized compilations of these samples have facilitated 
measurements of the relative expansion history across the 
redshift range 0 <z< 1 characterized by the equation-of-state 
parameter of dark energy (w = P/ (pc”)), and the measurement of 
the Hubble constant Hp, the current expansion rate determined 
from absolute distances. Measurements of w are constrained 
from the comparison of standardized SN Ia magnitudes over a 
wide range of redshifts obtained from different surveys with 
different observing-depth strategies. Measurements of Ho require 
very nearby (<50 Mpc, ~1 discovered per year) SNe Ia found 
by multiple surveys in galaxies that host calibrated primary 
distance indicators (e.g., Cepheids, tip of the red giant branch, 
TRGB), which are then compared to SNe in the Hubble flow, 
often from the same surveys. 

However, simply combining several subsamples into a large 
sample of SNe Ia does not provide meaningful gains without 
rigorous cross-calibration, self-consistent analysis of their light 
curves and redshifts, and characterization of their numerous 
sources of related uncertainties or covariance. As samples and 
compilations grow, ever greater attention must be paid to the 
control of systematic uncertainties, which would otherwise 
dominate sample uncertainties. 

This analysis, Pantheon+, is the successor to the original 
Pantheon analysis (Scolnic et al. 2018b) and builds on the 
analysis framework of the original Pantheon to combine an 
even larger number of SN Ia samples and include those that are 
in galaxies with measured Cepheid distances in order to be able 
to simultaneously constrain parameters describing the full 
expansion history (e.g., Qa, wo, and w,) with the local 
expansion rate (Hp). The original Pantheon compilation of 
1048 SNe Ia was used to measure a value (from SNe Ia alone) 
of w= — 1.090 + 0.220. Riess et al. (2016), in their measure- 
ment of the local expansion rate Ho, used a prerelease version 
of Pantheon based on Scolnic et al. (2015) and further 
augmented the sample as Pantheon did not extend to reach 
the low redshifts of the primary distance indicators at z < 0.01. 

Although there was significant overlap in data and analysis 
between the Pantheon measurement of w and the Riess et al. 
(2016) measurement of Hp, the Riess et al. (2016) measurement 
included several Cepheid-calibrator SNe Ia that were not 
included in Pantheon, and the fitting for Hp and parameters 
describing the expansion history were done independently 


rather than simultaneously. Dhawan et al. (2020) later 
established a framework for considering the covariance 
between SNe in primary distance indicator hosts and SNe in 
the Hubble flow. We build on that framework, which was 
developed originally for a redshift-binned Hubble diagram, and 
in this paper we create the first unbinned sample with 
covariance extending down to z= 0.001 that can be used to 
propagate correlated systematics for simultaneous measure- 
ments of Ho, Quy, Wo, and wy. We (i) analyze the largest set of 
cosmologically viable SN Ia light curves to date, (ii) include 
low-redshift samples to extend the lower bound in redshift to 
0.001, which contains the primary distance indicators (SNe in 
SHOES Cepheid host galaxies), (iii) propagate systematic 
uncertainties for both primary distance indicators and higher- 
redshift SNe simultaneously, and (iv) leverage the large strides 
made in the field of SNIa cosmology since the original 
Pantheon. 

This paper is the culmination of a series of papers that 
comprise the Pantheon+ analysis. A graphic of an overview of 
the numerous Pantheon+ supporting analyses, on which this 
paper heavily relies, is shown in Figure |. Details of each paper 
pertinent to this analysis are described in Section 3. One of 
these papers, Scolnic et al. (2022, hereafter S22), describes the 
sample of 1701 cosmologically viable SN Ia light curves of 
1550 distinct SNe, which we will refer to as “the Pantheon+ 
sample.” The redshifts and peculiar velocities of the SNe used 
here are given by Carr et al. (2021), and a comprehensive 
analysis of peculiar velocities is presented by Peterson et al. 
(2022). The cross-calibration of the different photometric 
systems used in this analysis can be found in Brout et al. (2022, 
hereafter Fragilistic), and calibration-related systematic uncer- 
tainty limits are determined by Brownsberger et al. (2021). The 
underlying SN Ia populations describing the data set are given 
by Popovic et al. (2021b). The model for intrinsic brightness 
variations was developed by Brout & Scolnic (2021) and then 
improved and evaluated by Popovic et al. (2021a). The novel 
systematic framework for simultaneous measurement of Ho and 
cosmology was developed by Dhawan et al. (2020), and an 
improved methodology for systematic uncertainties is 
described by Brout et al. (2021). 

In this work we discuss briefly the aforementioned papers in 
the context of their use in this analysis, evaluate several 
additional systematic uncertainties not addressed in these 
works, measure cosmological parameters, examine additional 
signals in the Hubble diagram, and compile systematic 
uncertainty budgets on cosmological parameters. A companion 
paper by the SHOES Team (Riess et al. 2022, hereafter R22) 
combines from this work 277 Hubble flow (0.023 < z< 0.15) 
SNe Ia and 42 SNe Ia in Cepheid-calibrator hosts, their relative 
distances, and their covariance, with the absolute distances of 
primary distance anchors (Cepheids and TRGB) from R22 in 
order to measure Hg under the assumption of flat ACDM. 
Similarly, in this work we utilize the full Pantheon+ sample of 
1550 SNela in combination with the R22 Cepheid host 
distances to show the impact of cosmological models with 
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Figure 1. Analysis roadmap of this work and supporting/complementary Pantheon+ and SHOES papers. Components of the analysis here are shown in blue. The 
companion paper (Riess et al. 2022, hereafter R22), which provides a constraint on Ho, requires the Hubble diagram and covariance computed in this work. Likewise, 
measurements of Ho in this work require the R22 Cepheid distance and covariance. Supporting papers are shown in gray boxes. 


more freedom than those used in R22 as well as the impact of 
SN-related systematic uncertainties on inferences of Hp. 

An important aspect of this work is the public release of the 
data and simulations used here that allow for the reproduction 
of multiple different stages of this analysis. In Appendix C, we 
present the numerous products that will be made available, 
including SN distances, redshifts, uncertainties, covariance, 
and extensive SNANA simulations (Kessler et al. 2009) of the 
data that model astrophysical effects, cosmological effects, and 
the observation/telescope effects of each survey down to the 
level of cadence, weather history, etc. We encourage the 
community to validate alternate analyses of the publicly 
released Pantheon+ sample on these simulations. 

The structure of the paper is as follows. In Section 2, we 
describe the methodology from fitting SN light curves to 
constraining cosmological parameters. Section 3 summarizes 
all of the inputs to the analysis including the data sample, 
calibration, and redshifts. In Section 4, we describe the 
cosmological results. Sections 5 and 6 are our discussions 
and conclusions, respectively. 


2. Methodology of Constraining Cosmological Parameters 
with SNe Ia 


2.1. Measuring Distances to SNe la 


To standardize the SN Ia brightnesses, we fit light curves 
using SNANA with the SALT2 model as originally developed 
by Guy et al. (2010) and updated in Brout et al. (2022), 
hereafter SALT2-B22. For each SN, the SALT2 light-curve fit 


returns four parameters: the light-curve amplitude x9 where 
mg = —2.5 log,9(xo); x1, the stretch parameter corresponding to 
light-curve width; c, the light-curve color, which includes 
contributions from both intrinsic color and dust; and fo, the time 
of peak brightness. Extinction due to Milky Way dust is 
accounted for in the SALT2 light-curve fitting. From the 
parameters mp, x;, and c, we standardize the SN brightnesses 
and infer distance moduli (/z), used in the Hubble diagram, with 
a modified version of the Tripp (1998) distance estimator. 
Following Kessler & Scolnic (2017, hereafter BBC), the 
distance modulus is defined as 


Spias Shosts (1) 


= mp + ax — Bc-— M 


where a and ( are global nuisance parameters relating stretch 
and color, respectively, to luminosity. M is the fiducial 
magnitude of an SN Ia, which can be calibrated by setting an 
absolute distance scale with primary distance anchors such as 
Cepheids. dp; is a correction term*> to account for selection 
biases that is determined from simulations following Popovic 
et al. (2021b), described in detail in Appendix A. djos¢ is the 
luminosity correction (step) for residual correlations between 
the standardized brightness of an SN Ia and the host-galaxy 


33 Dast analyses have the opposite sign + dpias; however, since the values of 
bias In the public release are meant to be subtracted, we change the sign 
compared to previous works. 
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where y is the magnitude of the SN Ia luminosity differences 
between SNe in high (M, > 10'° M5) and low (M, < 10'° M3) 
stellar-mass galaxies and where “hostless’ SNe have been 
assumed to reside in galaxies with low stellar mass. M, is the 
inferred stellar mass measured in units of solar mass (M5) from 
spectral energy distribution fitting to the photometry of each 
host galaxy, S is the step location (nominal analysis assumes 
S=10'°M.), and Ty, describes the width of the step. 

The total distance-modulus error, o,,, for SN i is described as 


2, 2 2 
Oi = Sf i, Cis M,,.i) O meas,i + O floor (Zis Cis M,.i) 
2 2! 2 
+ O'ens,i = O7i + Oo Vpec, i? (3) 


where Omeas is the measurement uncertainty of SALT2 light- 
curve fit parameters and their associated covariances (see 
Equation (3) of Kessler & Scolnic 2017) resulting from 
photometric uncertainties. The measurement uncertainty is 
scaled by f(z; c;, Mi) specific to each survey in order to 
account for selection effects that can reduce the observed 
scatter at the limits of each sample. The uncertainty contrib- 
ution from gravitational lensing as given by Joénsson et al. 
(2010) is Gens = 0.055 z. We note that, as discussed by Kessler 
et al. (2019a), the correct lensing distribution is utilized in 
simulations. The nominal distance-modulus uncertainty 
contribution due to the combination of redshift-measurement 
uncertainty (a,) and peculiar-velocity uncertainty (ype) have 
both been converted to distance-modulus uncertainty under the 
assumption of a cosmological model. Chen et al. (2022) noted 
that the optimal way to characterize redshift-measurement 
uncertainty at high redshifts (e.g., the Dark Energy Survey, 
DES, sample, z> 0.3) is to float the redshift and use the 
uncertainty in redshift as a prior in the light-curve fit. However, 
following previous analyses, we fix the redshift and include the 
associated distance-modulus uncertainty o, in Equation (3), 
which is a correct estimate at low redshifts (z < 0.1). Lastly, 
Ofloor represents the floor in standardizability owing to intrinsic 
unmodeled variations in SNe Ia such that 


2 2 2 
A floor (Zis Ci, M,.i) = OT scat (Zis Cis M,.i) +f O gray? (4) 


where o2.4:(zi, G;, My,;) is determined from a model that 


describes intrinsic brightness fluctuations, and ce is a single 


number representing a gray (color-independent) floor in 
standardizability for all SNe Ia; Bea is determined after the 
BBC fitting process in order to bring the Hubble diagram 
reduced xX? to unity. The details of o2..(z;, c;, M,.i), its model 
dependence, and its contribution to systematic uncertainties are 
discussed in further detail in Section 3.3.2 and Appendix A. 
To determine the distance-modulus values of all of the SNe, 
we follow the BBC fitting process with updates to increase the 
dimensionality of bias corrections in Popovic et al. (2021b). 
The likelihood (as given in Equation (6) of Kessler & 
Scolnic 2017) results in a cosmology-independent minimiza- 
tion of the free parameters (a, 3, y, and Ogray) that minimize the 
scatter in the Hubble diagram. While the BBC process was 
designed for utility for photometric cosmology analyses and 


Brout et al. 


uses SN Ia classification probabilities, the data analyzed here 
are a spectroscopically confirmed SN Ia sample, and therefore 
we set the non-Ia SN probabilities to zero for the whole sample. 


2.2. The Covariance Matrix 


Following Conley et al. (2011), we compute covariance 
matrices Ca, and Cyy, to account for statistical and systematic 
uncertainties and expected correlations between the SN Ia light 
curves in the sample when analyzing cosmological models. 
BBC produced both a redshift-binned and an unbinned Hubble 
diagram, enabling both binned and unbinned covariance 
matrices. For the original Pantheon (Scolnic et al. 2018b), 
Joint Light-Curve Analysis (JLA; Betoule et al. 2014), and 
DES3YR (Brout et al. 2019b), Cota and Cyyg, were redshift- 
binned matrices (or smoothed as a function of redshift) citing 
computational limitations. Following Brout et al. (2021), in this 
work we utilize the unbinned Hubble diagrams to create 
unbinned covariance matrices. The Pantheon+ sample (Scolnic 
et al. 2022) also includes “duplicate SNe Ia,” SNe Ia that have 
been observed simultaneously by numerous different surveys, 
so that statistical covariance Cyt, is computed as 


2 ae 
on b=J 


Cota @, J) = Caer + iene ’ (5) 
a + ee i ~ j and SN; = SN; 

where each row of the matrix corresponds to an SN light curve, 
the diagonal of Cstat is the full distance error (77) of the ith light 
curve, and where measurement noise from components other 
than the light curve itself is included as off-diagonal covariance 
between entries corresponding to light curves of the same SN 
(SN; = SN;) observed by two different surveys. 

Systematic uncertainties can manifest in three key places in 
the analysis: (1) from changing aspects affecting the light-curve 
fitting (e.g., survey photometry, calibration, SALT2 model), (2) 
from changing redshifts that propagate to changes in distance 
moduli relative to a cosmological model, and (3) from changes 
in the astrophysical or survey-dependent assumptions in the 
simulations used for bias corrections. For each of these 
categories, we examine all of the known significant sources 
of systematic uncertainty (7) with sizes S,, which result in 
residuals in the Hubble diagram relative to our baseline 
analysis (/gasg)- In order to compute the effect of systematics, 
we first define 


Ap, = iy = aes = (ret (Zs) ~~ Lovet (ZBASE))> (6) 


where Hy is the set of distances for systematic 7. For 
systematics that affect redshift, we have included a new 
methodology in Equation (6) that utilizes a reference 
cosmological model distance [;,¢(z) corresponding to flat 
ACDM (Qy = 0.3, Qa = 0.7). The f4,,5q(Z) and Ly,5q(ZBASE) 
are the cosmological model distances corresponding to red- 
shifts zy, and Zpasg. In order to propagate redshift effects into a 
distance x distance covariance matrix, the additional comp- 
ONneNt Lynog (Zi) — Umoa(ZBasE) accounts for the difference in 
inferred model distance. 

Assuming linearity between Aj,, and 7, we compute the 
derivative for each 7 in order to build a 1701 x 1701 
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which denotes the covariance between the ith and jth light- 
curve fit summed over the different sources of systematic 
uncertainty (~) with uncertainty o,, (see Section 3 for details). 
As shown by Brout et al. (2021), the o,, serve as priors on the 
known size of systematic uncertainties, but the data itself can 
constrain the impact of each systematic under the condition that 
information has not been collapsed by binning/smoothing (as 
was done for the original Pantheon, JLA, and DES3 YR). 

Fluctuations of the sample of light curves that pass the 
sample quality cuts (Table 2 of S22) for different systematics 
result in an ill-defined covariance matrix. To have a well- 
defined unbinned covariance matrix requires a subtle treatment 
in order to ensure that the sample is consistent in both the light- 
curve fitting and BBC stages across all systematics in the 
analysis. Quality cuts at the light-curve stage are only applied 
to the set of SNe based on their values found in the baseline 
analysis, and this SN sample is used for all systematic tests. We 
perform the BBC process twice—the first iteration to identify 
the subset of <1% of SNe for which bias corrections are unable 
to be computed, and a second iteration using only the common 
set of SNe that have valid bias corrections in all systematic 
variants. The final cosmology sample of 1701 light curves that 
satisfy all criteria is described in detail in S22 (see the 
“Systematics” row in Table 2 of $22). 

Finally, the statistical and systematic covariance matrices are 
combined and used to constrain cosmological models: 


Cytat-+syst = Cotat =F Cyst: (8) 


2.3. Cosmology 


Constraining cosmological models with SN data using y* 
has been used in previous SN Ia cosmology analyses (e.g., 
Riess et al. 1998; Astier et al. 2006) and first included 
systematic covariance in Conley et al. (2011). Here we follow 
closely the formalism of Conley et al. (2011) where 
cosmological parameters are constrained by minimizing a a 
likelihood: 


—21n(L) => x? => ADT Crit AD, (9) 


where D is the vector of 1701 SN distance-modulus residuals 
computed as 


AD; = Hi — Lmodei (Zi): (10) 


and each SN distance (j1;) is compared to the predicted model 
distance given the measured SN /host redshift ((4modei(Z;)). The 
model distances are defined as 


Hmodei (Zi) = 5 log(dr(zi)/10 pe), (11) 
where d; is the model-based luminosity distance that includes 
the parameters describing the expansion history H(z). For a flat 
cosmology (QQ; = 0), the luminosity distance is described by 

dz! 
H(z!) 
where d;(z) is calculated at each step of the cosmological fitting 
process, and the parameterization of the expansion history 


dy(z) = (1 + ae f (12) 
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(used in Equation (12) and therefore in the likelihood 
Equation (9)) in this work is defined as 


H(z) = Hy VQu0 + 32+ KAA + Z3G™, (13) 


See Hogg (1999) for the forms of the expansion history H(z) 
used in the case that the assumption of flatness is relaxed. 

The parameters M (Equation (1)) and Ho (Equation (13)) are 
degenerate when analyzing SNe alone. However, we also 
present constraints that include the recently released SHOES 
Cepheid host distance anchors (R22) in the likelihood that 
facilitates constraints on both M and Hp. 

When utilizing SHOES Cepheid host distances, the SN 
distance residuals are modified as follows: 


Ap!. = Ll; — pee 
~~ . 
Li — Umode (Zi) otherwise, 


i € Cepheid hosts (14) 


Cepheid 


where [U; is the Cepheid calibrated host-galaxy distance 


provided by SHOES, and where pu; — poe is sensitive to the 
parameters M and Hp and is largely insensitive to Qy, or w. We 
also include the SHOES Cepheid host-distance covariance 
matrix Ca) presented by R22 such that the likelihood 
becomes 

—2In(L’) = AD! (Cyitssyst + Cotatrsyst) | AD', (15) 
where ce +syst denotes the SN covariance. 

We evaluate the likelihoods with the PolyChord (Handley 
et al. 2015) sampler in the CosmoSIS package (Zuntz et al. 
2015) using 250 live points, 30 repeats, and an evidence 
tolerance requirement of 0.1. This resulted in converged chains 
containing 1000-3000 independent samples. We verified the 
SN-only results with CosmoMC (Lewis & Bridle 2002) and 
with the fast cosmology grid-search program in SNANA. The 
likelihood for Pantheon+ and R22 Cepheid host distance 
samples will be made available in the public version of 
CosmoSIS. In this work we also utilize the additional public 
likelihoods in CosmoSIS in order to combine with and assess 
agreement with external cosmological probes: Planck (Planck 
Collaboration et al. 2020) and baryon acoustic oscillations 
(BAOs; likelihoods discussed in Section 4). 


In this work we _ investigate four cosmological 
parameterizations: 
1. Flat ACDM: Quy, is floated, and we fix w=—1 and 
Qu a QA —— 1. 


2. ACDM: Qy and Q, are floated, and we fix w=—1. 

3. Flat wCDM: w and Qy are floated, and we fix 
QytQya=1. 

4. Flat wow,CDM: w=wo+w,1+z), Qu, Wo, Wa are 
floated, and we fix Qy+Q,=1. 


We blind our analysis in two ways simultaneously. First, we 
blind the binned distance residuals output by the BBC fit, as 
cosmological parameters could be inferred visually from 
simply looking at the Hubble diagram. Second, in order to 
prevent accidental viewing of the cosmological parameters 
themselves, the CosmoSIS chains were shifted by unknown 
values following the formalism of Hinton (2016). 


3. Data and Analysis Inputs 


Here we review each component of the data set and analysis. 
We discuss the fundamental purpose, the baseline treatment in 
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this analysis, and the systematic uncertainties associated with 
each aspect (if applicable). The impact of systematics in both 
distance and cosmological inference is shown in Section 4. We 
provide a brief overview of this section here. 

Data 


Section 3.1.1: SN Ia Light Curves 
Section 3.1.2: Redshifts 

Section 3.1.3: Peculiar Velocities 
Section 3.1.4: Host-galaxy Properties 


Calibration and Light-curve Fitting 


Section 3.2.1: Calibration 
Section 3.2.2: SALT2 Model 
Section 3.2.3: Milky Way Extinction 


Simulations 


Section 3.3.1: Survey Modeling 
Section 3.3.2: Intrinsic Scatter Models 
Section 3.3.3: Uncertainty Modeling 
Section 3.3.4: Validation 


3.1. Data 
3.1.1. SN la Light Curves 


Purpose: The flux-calibrated light-curve photometry is fit to 
determine the SALT2 parameters used in standardization 
(Equation (1)). 

Baseline: The light-curve data is described in detail by S22 
and references therein. The full set of spectroscopically 
classified photometric light curves is compiled from 18 
different publicly available and privately released samples. In 
total, 2077 SN light-curve fits converged using SALT2; after 
quality cuts are applied (Table 2 of $22), this results in 1701 
SN light curves of 1550 unique SNe Ia usable for cosmological 
constraints. The sample includes a 3.50 Hubble residual outlier 
cut to remove five potential contaminants that are likely 
nonnormal Type Ja or misidentified redshifts. The sample of 
cosmologically viable light curves includes 81 light curves of 
42 SNe used to calibrate Cepheid brightnesses as utilized 
by R22. The survey SN photometry compiled in Scolnic et al. 
(2022) and analyzed here is from DES* (Brout et al. 2019a; 
Smith et al. 2020a), Foundation! (Foley et al. 2018), Pan- 
STARRS (PS1; Scolnic et al. 2018b), Supernova Legacy 
Survey (SNLS; Betoule et al. 2014), Sloan Digital Sky Survey 
(SDSS; Sako et al. 2011), Hubble Space Telescope (HST; 
Gilliland et al. 1999; Riess et al. 2001, 2004, 2007; Suzuki 
et al. 2012; Riess et al. 2018), and Low-z (grouped together as 
LOSS_1'; Ganeshalingam et al. 2010; LOSS_2'; Stahl et al. 
2019; SOUSA’*?; Brown et al. 2014; CNIa0.02'; Chen et al. 
2020; CSP; Krisciunas et al. 2017b; CfA1; Riess et al. 1999; 
CfA2; Jha et al. 2006; CfA3; Hicken et al. 2009; CfA4; Hicken 
et al. 2012; and numerous smaller low-redshift samples! of one 
to two SNe given by Milne et al. 2010; Stritzinger et al. 2010; 
Tsvetkov & Elenin 2010; Zhang et al. 2010; Krisciunas et al. 
2017a; Burns et al. 2018; Gall et al. 2018; Burns et al. 2020; 
Kawabata et al. 2020.) 

Systematics: See Calibration (Section 3.2.1). 


34 Not included in Pantheon 2018. 
39 https: //pbrown801.github.io/SOUSA/ 
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3.1.2. Redshifts 


Purpose: The peculiar-velocity corrected cosmic microwave 
background (CMB) frame redshift of each SN/host is required 
to compare the inferred distance to a distance predicted by a 
cosmological model, as given in Equation (10). Additionally, 
heliocentric redshifts are required in the SALT2 light-curve fits 
in order to shift the model spectrum to match the data. 

Baseline: The redshifts for all of the SNe (and their host 
galaxies, depending on what is available) are provided by Carr 
et al. (2021), who performed a comprehensive review of 
redshifts for the Pantheon+ samples and made numerous 
corrections. Carr et al. (2021) reported the heliocentric redshifts 
for each SN and converted the redshift into the CMB frame. 
The redshifts of the Pantheon+ sample cover a range of 
0.001 <z<2.3. While redshifts of the 42 Cepheid host 
calibrator SNe are included, they are not used in the 
comparison of SN Ia magnitudes to the Cepheid distance scale 
and are only provided for reference and for SALT2 fitting. 

Systematics: Following Carr et al. (2021), we apply a 
coherent shift to each redshift of +4 x 10-°. This was 
conservatively stated by Calcino & Davis (2017) for the 
potential size of a local void bias and by Davis et al. (2019) as a 
potential measurement bias. 


3.1.3. Peculiar Velocities 


Purpose: Peculiar motions of galaxies arise from coherent 
flows, motion of halos, inflow into clusters or superclusters, 
and intragroup motion. Corrections are applied to the observed 
redshifts (after light-curve fitting) based on peculiar-velocity 
maps derived from independent large spectroscopic galaxy 
surveys. 

Baseline: The nominal peculiar velocities used for this 
analysis were determined by Peterson et al. (2022) from a 
comparison of multiple treatments of peculiar-velocity maps and 
group catalogs. Corrections were applied by Carr et al. (2021) 
for the Pantheon+ sample. The baseline corrections are based on 
2M+-+ (Carrick et al. 2015) with gl obal parameters found in 
Said et al. (2020) and combined with group velocities estimated 
from Tully (2015) group assignments. The oypec in Equation (3) 
is found using 240kms~' after accounting for uncertainties 
propagated into the covariance matrix described below. This 
Ovpec floor is in agreement with what was used in Peterson et al. 
(2022), and for the SNe between 0.001 < z< 0.02, it is likely a 
conservative estimate, as Kenworthy et al. (2022) found a floor 
of 155+25kms_' for the most nearby SN calibrators. This 
apparent reduction at the lowest redshifts may be due to the 
peculiar-velocity maps having higher fidelity at these redshifts 
and because Pantheon+ has relatively better virial-group 
information at these redshifts. 

Systematics: Peterson et al. (2022) discuss multiple viable 
alternatives for the treatment of peculiar velocities. The first 
approach is to use the 2M+-+ corrections (Carrick et al. 2015) 
integrating over the line-of-sight relation (LOS) between distance 
and the measured redshift. We take this variation as the first 
systematic with 0%, = 0.5. The second approach is to use 
the Two-Micron All-Sky Redshift Survey (2MRS; Lilow & 
Nusser 2021) peculiar-velocity map; however, differences 
between 2MRS and 2M+-+ at very low redshift (z < 0.01) cause 
numerical stability issues for off-diagonal C,, elements. We 
incorporate only the diagonal differences between 2MRS and 2M 
++ into Cyys with oF = 0.5. As a numerically stable estimate of 
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the 

off-diagonal terms, we use the 2M+-+ velocities transformed by 
the slope and offset difference between the 2M+-+ and 2MRS 
maps found in Peterson et al. (2022). The two approaches added 
in quadrature result in an effective a4, = 1.0. 


3.1.4. Host-galaxy Properties 


Purpose: The observed host-galaxy mass versus SN luminosity 
relation is used to standardize the SN Ia brightnesses in two ways. 
First, simulations of the data set include correlations between SN 
color and SN stretch and host properties such as dust as a function 
of host mass, following Popovic et al. (2021a). Second, a further 
residual correction is applied in the Tripp Equation (1) where the 
“mass step” y is fit in the BBC stage. 

Baseline: The host-galaxy stellar masses are presented by $22 
and references therein. Masses are determined for all host 
galaxies, and star formation rates and morphologies are also 
included the low-z sample. In the baseline analysis, we apply the 
mass step at 10!° M., following Pantheon and DES3YR. 

Systematics: Several independent analyses (Sullivan et al. 
2010; Childress et al. 2013a; Kelsey et al. 2021) have 
suggested that the optimal location of the mass step could 
range between 10° Mz and ig M>. We therefore include a 
systematic uncertainty where the mass step occurs at io" Mo. 


3.2. Calibration and Light-curve Fitting 
3.2.1. Calibration 


Purpose: Photometric calibration of each passband in each 
survey is needed to fit light curves and facilitate comparison of 
the brightnesses of SNe across different telescopes /instru- 
ments/filters. Photometric calibration is also important to 
homogenize spectrophotometric data sets used in the SALT2 
model training. 

Baseline: The calibration of all 25 photometric systems used 
in this work is discussed in Brout et al. (2022). The outputs of 
Fragilistic are a best-fit calibration solution for each of the 105 
passbands and a joint 105 x 105 covariance matrix that 
describes the covariance between the zero-point calibrations 
of all passbands that arise from using a single common stellar 
catalog to tie all surveys together (PS1). 

Systematics: The systematics due to calibration and their 
impact are discussed in detail in Fragilistic. We estimate the 
impact of the correlated filter zero-point and central wavelength 
uncertainties by refitting SALT2 light curves (with retrained 
SALT2 models; see Section 3.2.2) using nine realizations of 
the 105 zero-points. For each of the nine realizations, a value of 
o%, =1 i 9 is adopted such that they add in quadrature to ~1. 
The uncertainty in modeling the spectrum of the HST primary 
standard star C26202 has been tripled to account for the recent 
update in Bohlin et al. (2020); it is now set to 15 mmag over 
7000 A (oy,=3 for a systematic of 5mmag over 7000 A). 
Lastly, an additional conservative systematic is included only 
for the CSP SNe to account for the 2% recalibration in CSP 
tertiary stellar magnitudes from Stritzinger et al. (2010) to 
Krisciunas et al. (2017b; oy, = 1). 


3.2.2. SALT2 Model 


Purpose: The trained SALT2 model is required to fit light 
curves and determine the light-curve parameters (m,, c, and x;) 
for each SN used in Equation (1). 
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Baseline: We use the Fragilistic calibration solution and 
newly trained SALT2-B22 model,*° which was developed 
following the formalism of Guy et al. (2010) and Taylor et al. 
(2021). The SALT2 model includes a component of training 
statistical uncertainty, which is incorporated in the fitted light- 
curve parameters. 

Systematics: For each of the nine correlated realizations of 
Fragilistic filter zero-points and central wavelengths discussed 
above (for calibration), we simultaneously retrain the SALT2 
model. Additionally, to conservatively account for a possible 
systematic from the redevelopment of the SALT2 model- 
training process itself, we adopt an additional systematic by 
fitting the data set with the SALT2 model trained by Betoule 
et al. (2014) and applying a scaling of 7, = 1/3 (see Section 5 
and Figure 15 for impact). 


3.2.3. Milky Way Extinction 


Purpose: Values of the Milky Way (MW) dust extinction, 
E(B — V)mw, are applied to the SALT2 model spectra during 
both the model-training process and during the data light-curve 
fitting process. The “extinction curve” describes the relation 
between the amount of reddening and extinction as a function 
of wavelength. 

Baseline: We account for MW extinction using maps from 
Schlegel et al. (1998), with a scale of 0.86 following Schlafly 
et al. (2010). We assume the MW extinction curve from 
Fitzpatrick (1999) with Ry = 3.1. 

Systematics: Similarly to Pantheon, we adopt a global 5% 
uncertainty scaling of E(B—V)yw based on the fact that 
Schlafly & Finkbeiner (2011), in a reanalysis of Schlafly et al. 
(2010), derive smaller values of reddening by 4%, despite 
using a very similar SDSS footprint (a, = 1). While Schlafly & 
Finkbeiner (2011) found that their results prefer the Fitzpatrick 
(1999) extinction curve, we conservatively include an addi- 
tional systematic uncertainty in the MW extinction curve and 
analyze the data (training and light-curve fit) using the 
extinction curve from Cardelli et al. (1989) and apply a 
systematic scaling of o,, = 1/3, as this reflects the preference of 
Fitzpatrick (1999) over Cardelli et al. (1989). 


3.3. Simulations 
3.3.1. Survey Modeling 


Purpose: We utilize catalog-level simulations of large 
samples of SN Ia (> 1,000,000 per survey) light curves. SNANA 
simulations specific to each survey in our analysis are 
prescribed by each aspect of acquiring an SN Ia sample. As 
detailed in Figure 1 of Kessler et al. (2019a), the simulations 
require three main sets of inputs: 

A source model for generating SNe with realistic astrophysical 
properties and applying cosmological effects such as redshifting, 
dimming, lensing, peculiar velocities, and MW extinction. 

A noise model, unique to each survey, for applying 
instrumental and atmospheric noise to determine a detection 
efficiency (““DETEFF’). 

A trigger model, unique to each survey, that includes the 
observing cadence and describes an efficiency as a function of 
B-band peak magnitude for detecting SNe and obtaining a 
spectroscopic confirmation (“SPECEFFP’). 


36 Released publicly at pantheonplusshOes.github.io. 
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Figure 2. Comparison between observed data (black points) and simulations 

(blue lines) for the largest subsamples in this analysis: DES, HST, SDSS, 

SNLS, PS1, LOW-z, and Foundation (FOUND). We compare three key 

distributions: the SALT2 light-curve-fit parameters x, and c are shown as well 

as the measured redshift. 


These simulations for each survey are combined and used to 
forward model the underlying populations of the SN properties 
(Popovic et al. 2021a, 2021b) and to determine the expected biases 
in measured SN distances that follow from the known selection 
effects. These biases are corrected in the 6,;,, term of Equation (1). 

Baseline: Depicted in Figure 2 are the distributions of the 
key observables (z, x,, and c) for both data and simulations of 
each survey used in this analysis. We find good agreement 
between the data and simulations, as described in detail by 
Popovic et al. (2021a) and Popovic et al. (2021b). We note that 
the agreement in the redshift dimension is achieved despite not 
explicitly tuning the redshift distribution of surveys. 

We simulate SNe in LOW-z and Foundation down to z = 0.001. 
Novel for this work specifically are the simulations of primary 
distance indicator hosts of SNe in the range 0.001 <z< 0.01, 
which are assumed to have the same color and stretch populations 
as those of their respective surveys (LOW-z and FOUND), and 
specifically over this redshift range, they are assumed to be 
complete with flat spectroscopic selection efficiency. These 
simulations facilitate bias corrections to the Cepheid-calibrator 
SNe and thus the propagation of modeling systematics to the SNe 
used in the companion SHOES analysis (R22). 
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The simulation inputs for survey cadence, DETEFF, and 
SPECEFF functions have been evaluated in many analyses 
over the past decade. Table 1 shows a summary of where we 
obtain these inputs for each survey. Survey metadata is used to 
model the cadence and instrumental properties, if available, 
such as for FOUND, SDSS, PS1, DES, and SNLS. LOW-z data 
do not provide such metadata, and thus the cadence and noise 
properties are extracted from the data as described in Kessler 
et al. (2019a) following the procedure developed by Scolnic 
et al. (2018b), which assumes that the LOW-z subset of SNe is 
magnitude-limited. These are simulations of the CfA and CSP 
samples, but not of the newer samples included in this work 
(LOSS, SOUSA, and CNIa0.02), thereby implicitly assuming 
that the CfA and CSP samples have similar selection effects 
and therefore distance biases as the newer additions. To 
simulate SN-host correlations, a catalog of host-galaxy proper- 
ties and specifically their stellar-mass distributions is taken 
from Popovic et al. (2021b). The simulations used for bias 
corrections for all surveys are performed in ACDM (w = — 1.0, 
Qy = 0.3, XQ, =0.7) with the SALT2-B22 model. 

Systematics: We increase the S/N of each simulation by 
20%, resulting in all survey simulated distributions changing 
by more than lo, as a single conservative systematic in the 
determination of the selection biases. Kessler & Scolnic (2017) 
showed that the sensitivity of the bias corrections to the input 
cosmology is relatively weak; this was confirmed by Brout 
et al. (2019b) and found to be a negligible contribution to SN Ia 
uncertainty budgets. We therefore do not include this as an 
additional systematic. 


3.3.2. Intrinsic Scatter Models 


Purpose: A model of the intrinsic SN brightness variations, 
called “intrinsic scatter,” is needed to account for the observed 
residual variation in SN Ia standardized luminosities that 
exceeds expectations from measurement uncertainties alone. 
In addition, models of the true (“parent”) populations of SN Ia 
SALT2 parameters c and x; are required for the source model 
in SNANA. The intrinsic scatter model is utilized in the bias- 
correction simulations. 

Baseline: We utilize the Brout & Scolnic (2021, 
hereafter BS21) model that prescribes SN Ia scatter into two 
color-dependent components: (i) a standard cosmological color 
law specific to SNe Ia and (ii) additional dust-based color laws 
and dust extinctions that vary with each galaxy/SN. This 
approach is preferred because of its novel replication of the 
observed relationships between SN color and residual Hubble 
diagram scatter as well as its ability to replicate the “mass step” 
as a function of SNIa color. We use the scatter model 
parameters from BS21 with improvements from Popovic et al. 
(2021a) in our baseline bias-correction simulations; because 
the BS21 model includes within it the parent c population, we 
also utilize the separate parent population for x, derived by 
Popovic et al. (2021b). Improving upon Scolnic & Kessler 
(2016), Popovic et al. (2021b) fit for parent populations in bins 
of mass to account for host—SN Ia relationships. Popovic et al. 
(2021b) split their populations into high- and low-redshift 
groups, and notably for low-redshift surveys, the x; populations 
are fitted with a two-Gaussian model to recreate the observed 
double peak in the x, distribution. 

Systematics: We include two categories of systematics for 
the intrinsic scatter model and parent populations: (1) 
different models of intrinsic scatter, and (2) determination 
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Table 1 
References for Inputs to SNANA Simulations Used for This Analysis 
Survey Cadence DETEFF SPECEFF 
LOW-z Scolnic et al. (2018b) Kessler et al. (2019b) Scolnic et al. (2018b) 
FOUND Jones et al. (2019) N/A Jones et al. (2019) 
SDSS Kessler et al. (2013) Kessler et al. (2009) Popovic et al. (2021b) 
PS1 Jones et al. (2018a) Jones et al. (2018b) Scolnic et al. (2018b) 
DES Smith et al. (2020b) Kessler et al. (2015) Abbott et al. (2019) 
SNLS Kessler et al. (2013) N/A Popovic et al. (2021b) 
N/A N/A 


HST Scolnic et al. (2018b) 


Note. We give references for the “Cadence,” which describes the observing history; the “DETEFF,” which describes the detection efficiency based on the signal-to- 
noise ratio (S/N); and the “SPECEFF,” which describes the spectroscopic selection efficiency as a function of SN magnitude. 


of parameters for the BS21 model. For the former, we use 
two additional scatter models from Kessler et al. (2013) that 
have been used in previous cosmology analyses (JLA, 
Pantheon, and DES3YR). These are (1) the “G10” model 
based on Guy et al. (2010), which describes ~70% of the 
excess Hubble scatter from “gray” variations and the 
remaining scatter from wavelength-dependent variations, 
and (2) the “C11” model based on Chotard (2011), which 
describes ~30% of the excess Hubble scatter from coherent 
variations, and the remaining scatter from wavelength- 
dependent variations. For the G10 and C11 scatter models, 
bias corrections are performed in seven dimensions as given 
by Popovic et al. (2021b). For the systematic uncertainty in 
the determination of the BS21 model parameters, we adopt 
three different viable sets of dust and intrinsic SN popula- 
tions from Popovic et al. (2021a). These populations are the 
best-fit (maximum likelihood) parameters (hereafter P21), the 
mean posterior set of parameters, and a set that represents a 
lo fluctuation in the uncertainty. Lastly, while the BS21 and 
P21 models impact the simulated bias corrections, the SALT2 
training and light-curve fitting has not been altered. The 
choice of scatter model is propagated through the simulations 
used for the bias corrections applied in Equation (1) and for 
the uncertainty modeling in o,c¢a, of Equation (4). 


3.3.3. Distance-modulus Uncertainty Modeling 


Purpose: To match the reported SN _ distance-modulus 
uncertainties (Equation (3)) to the scatter in distance that is 
observed in the data. 

Baseline: The BS21 model parameters have been fit to the 
observed scatter in the data set. We can utilize large BS21 
simulations to determine Oycat(z, c, M,) after accounting for 
selection effects. The efficacy of this method is shown in 
Figure 3, which demonstrates good agreement between the 
observed rms of the Hubble residuals and the uncertainties of 
the distance-modulus values. 

Systematics: To conservatively account for how SN cosmol- 
ogy was done in the past (JLA and Pantheon), in Equation (3) 
we set Oscat(z, c, M,)=0 and allow only a single Ogray 
parameter to replicate the methodology used with historic 
intrinsic scatter models (G10 and C11). However, we note that 
for G10 and C11, the trends in rms seen for the data in Figure 3 
do not match the reported uncertainties. 


3.3.4. Validation 


Purpose: To verify that our analysis can recover input 
values in data-sized simulated samples and does not produce 


— log(M,/Mza)> 10: RMS zu 
— log(M,/Ma)< 10: RMS pu 
---- log(M,/Mz)> 10: Mean o, 
---- log(M,/Mz)< 10: Mean o, 
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Figure 3. Pantheon+ distance-modulus uncertainties (shown as the dashed 
lines with mean g,, and split on host mass) in comparison to the observed rms 
of the distance-modulus residuals (shown as solid lines as rms ju split on host 
mass), as a function of color. This shows that the distance errors are adequately 
modeled (Equation (4)) as a function of SN color and host stellar mass. In 
previous analyses, the uncertainties were roughly flat as a function of color. 


biases. Such tests are sensitive to the light-curve fitting and 
BBC technique (as well as implementation and coding 
errors); however, they are not sensitive to certain aspects of 
the analysis such as the assumption of the SALT2 model or 
photometric calibration. 

Baseline: We perform an end-to-end test of our baseline 
analysis pipeline from survey photometry catalog-level 
simulations. We create 20 realizations of each survey in an 
arbitrary cosmological model (w= -—1): 10 with the BS21 
scatter model and 10 with the G10 scatter model. We perform 
light-curve fitting, apply bias corrections, compile into 10 
Hubble diagrams, and maximize the cosmological likelihoods 
(Equation (9)) using a fast cosmology grid-search program in 
SNANA (Kessler et al. 2009), with approximate priors from 
CMB measurements (Planck Collaboration et al. 2020) to 
obtain best-fit cosmological parameters and uncertainties. For 
the BS21 model simulations, we recover a mean best-fit 
w = —1.012 + 0.011, and for the G10 model simulations, we 
recover a mean best-fit w = —0.983 + 0.015; both are within 
~lo of the input cosmology. The 20 realizations are made 
available publicly*’ along with bias-correction simulations. 


37 Will be made available after publication at pantheonplusshOes.github.io. 
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Table 2 
Standardization Parameters and Results 


BBC Fit CosmoSIS Fit 
Model a B Ogray Y rms In(L) 
BS21 0.148(4) 3.09(4) 0.00 —0.003(7) 0.171 —1635 
P21 0.145(5) 3.00(5) 0.00 0.019(10) 0.171 —1674 
G10 0.153(4) 2.98(5) 0.10 0.054(7) 0.173 —1676 
Ci 0.153(4) 3.44(6) 0.12 0.053(8) 0.173 —1681 


Note. The nuisance parameters, as defined in Equations (1) and (3) are given 
here for different assumptions about the intrinsic scatter model, as described in 
Section 3 (intrinsic scatter model). The fact that cin ~0O and y~0 for 
the BS21 and P21 models is due to modeling the scatter and mass step as part 
of the BBC process, which is discussed in further detail in Appendix A. 
The BS21 is the baseline choice for intrinsic scatter. The rms is given in units 
of magnitudes. The Hubble diagram likelihood values for each model (L) 
include an uncertainty normalization term. 


4, Results 
4.1. Standardization Parameters 


The standardization nuisance parameters a, 3, 7, and Ogray 
defined in Equations (1) and (3) are shown in Table 2 for each 
of the scatter models used in this work. The best-fit a are 
similar across scatter models to within ~ lo. The best-fit @ 
values differ across models owing to different treatments of 
SN Ia color; however, the values for the baseline dust model 
(BS21) and the P21 dust model are self-consistent. 

As shown in Table 2, the additional @g;ay term for the BS21 and 
P21 models is found to be zero. As discussed in Section 2, this is 
consistent with the expectation that if the simulations correctly 
model the intrinsic scatter and noise of the data, the Ogca(z, c, M,) 
term of Equation (3) is sufficient to describe the distance-modulus 
uncertainties with gray =0. As discussed in Appendix A, for 
our G10 and C11 systematic treatment, o.ca(z, c, M,) is set to 0, 
and therefore Og;ay * 0.10 approximates the scatter, though it does 
not account for the observed color dependence. 

Table 2 also shows that the best-fit host stellar-mass corrections 
(y) are consistent with zero for BS21 and P21. This is in 
agreement with the findings of Popovic et al. (2021a), that 
modeling the intrinsic scatter in bias-correction simulations with 
correlations that match those in the observed data removes the 
need for ad hoc corrections in intrinsic brightness (i.e., y= 0). 
This can also be seen in Figure 5. For the bias correction based on 
the GIO and Cll models that do not include any mass 
dependence, the resulting y is ~0.05 found at 7o confidence. 


4.2. The Hubble Diagram and Distance Covariance Matrix 
4.2.1. The Hubble Diagram 


The Pantheon+ Hubble diagram of 1701 SN Ia light curves 
compiled from 18 different surveys and ranging in redshift 
from 0.001—2.26 is shown in the top panel of Figure 4. The 
bottom panel of Figure 4 shows the residuals to the best-fit 
cosmology (Equation (10)). Best-fit cosmological parameters 
will be presented in the following subsections. 

Shown in Table 2 is the total observed scatter (rms) in the 
Hubble diagram residuals to the best-fit model (bottom of Figure 4) 
for different scatter models. The BS21 model results in the lowest 
Hubble diagram rms and x, a>5o improvement determined 
from the difference in likelihoods relative to the G10 and Cll 
scatter models. Additionally, BS21 results in no significant 
residuals in the Hubble diagram as a function of SALT2 c, SALT2 
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x;, and host galaxy mass (Figure 5). The observed scatter of 
~0.17 mag is larger than that seen in the original Pantheon, 
because Pantheon+ extends to lower redshifts and thus is more 
impacted by scatter induced by peculiar velocities. If we set the 
minimum redshift to 0.01, the total scatter is reduced to 0.15 mag, 
matching that of Pantheon. Finally, compared to the original BS21 
analysis, P21 uses a more rigorous fitting process that is optimized 
to better characterize SN Ia colors and intrinsic scatter in addition 
to Hubble residuals. For this reason, the improvements of P21 are 
not solely described by the cosmological model likelihood £ of 
Table 2. We therefore have included the use of P21 population 
parameters as a systematic uncertainty. 


4.2.2. The Very Nearby Hubble Diagram 


We note from Figure 4 that in the very nearby universe, 
zZ< 0.008 (v < 2400 km ), the mean of the Hubble diagram 
residuals is positive by ~5% at ~ 20 significance. This is seen 
after the use of peculiar-velocity maps from either 2M+-+ or 
2MRS. A similar signal is also seen in the Hubble residuals of 
the Cepheid distances (Kenworthy et al. 2022). A bias of 
roughly this size and direction is expected in the presence of 
measurement errors and unmodeled peculiar velocities, which 
scatter more objects down from higher redshifts and a greater 
volume than from the reverse. This effect is significant only for 
the most nearby galaxies (z < 0.008). In Figure 4, we include 
the prediction (dashed line) for this bias assuming 250kms_' 
uncorrected velocity scatter (not a fit). 

In the three-rung distance ladder utilized to measure Hp by the 
SHOES Team (R22) and in Equation (14) in this work, the nearby 
(z < ~0.01) Hubble diagram is not used. Rather, only the distance 
moduli from such nearby SNe are used in the SN-Cepheid 
absolute distance calibration in the second rung. Furthermore, in 
the R22 measurement of the Hubble flow, only SNe with redshifts 
z > 0.023 are used in the third rung to limit sensitivity to peculiar 
velocities. This approach is insensitive to the volumetric redshift 
scatter effects, and there is no resulting impact on the R22 Ab. 
However, more local measurements of Hy from, for example, a 
two-rung distance ladder using primary distance indicators like 
Cepheids and TRGB and their host redshifts (mostly at z < 0.01), 
are more sensitive to peculiar velocities and the volumetric bias 
they induce, and are likely to be biased low at the few percent 
level if not appropriately accounting for this expected bias 
(Kenworthy et al. 2022). For measurements of other cosmological 
parameters (e.g., w or Qyy) with Pantheon+ described in the 
following subsections, the mean Hubble residual biases of the 
Low-z and Foundation sample are ~2mmag and ~1 mmag, 
respectively, and are considered to be negligible. 


4.2.3. The Distance Covariance Matrix 


Built following Equation (7), the 1701 x 1701 systematic 
distance covariance matrix is shown in Figure 6. The sample is 
sorted by survey and redshift to help visualize the covariances. 
The Hubble diagram residuals (Equation (10)) that are used to 
build the covariance matrix are shown in Figure 7 for several 
example sources of systematic uncertainty. As discussed in 
Appendix C, the information used to create the Hubble diagram as 
well as the covariance matrix is publicly available,** and tools to 
read in this information are in CosmoSIS. The SDSS 
subsample contributions to the covariance matrix (Figure 6) 


oe pantheonplusshOes.github.io 
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Figure 4. Top panel: the Pantheon+ “Hubble diagram” showing the distance modulus ju vs. redshift z. The 18 different surveys are each given different colors. Bottom 
panel: the distance-modulus residuals relative to a best-fit cosmological model with binned data for reference (black points). Both the data errors and the binned data 


errors include only statistical uncertainties. At z < 0.01, the sensitivity of peculiar velocities is very large, and the uncertainties shown reflect this uncertainty. The 
dashed line is the predicted Hubble residual bias stemming from biased redshifts due to volumetric effects in the very nearby universe (assuming 250 km s 


uncorrected velocity scatter). 


* — Step at Median: -0.002 + 0.007 * — Step at Median: 0,006 + 0.007 - ¢ — Step at Median: -0.002 + 0,007 


LL — Lmodel 


0.0 0. 


SALT2 c 


—0.2 


Figure 5. Pantheon+ sample Hubble diagram residuals (teal) to the best-fit 
cosmology (4 — Umode1) for the baseline analysis as a function of SALT2 c, 
SALT2 x;, and host-galaxy stellar mass M,. Distances (4) follow Equation (1) 
and include a, 2, Spias, aNd Sos: Corrections. Binned data are shown for 
reference (black). No significant residual correlations are seen. 


stand out visually due to their strong spectroscopic selection 
function. 


4.3. Constraints on Cosmological Parameters from Pantheon+ 
and SHOES 


Parameter constraints from the Pantheon+ SNelIa and 
SHOES Cepheid host absolute distances are shown in Table 3 
for flat ACDM, ACDM, flat wCDM, and flat wow,CDM. 
Unless otherwise stated, constraints on cosmological para- 
meters include both statistical and systematic uncertainties. 
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Figure 6. The systematic covariance matrix as defined in Equation (7). To 
show the inherent structure, the data set is sorted by survey and within each 
survey (colored boxes), by redshift. “CALIB” are the set of 81 SN light curves 
in the SHOES Cepheid-calibrator galaxies. The shading corresponds to the size 
of the covariance in magnitudes. 


From the Pantheon+ SNe Ia, for a flat ACDM model, we find 
Ou = 0.33440.018. We note that SHOES (R22) utilizes 
Pantheon+ SNe at z < 0.8 to constrain the deceleration parameter 
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Table 3 
Results for Cosmological Models 
Qu On Ho Wo Wa 
Pantheon+ and SHOES: All Models 
Flat ACDM 0.334 + 0.018 0.666 + 0.018 73.6+ 1.1 
ACDM 0.306 + 0.057 0.625 + 0.084 T4411 vee 
Flat wCDM 0.30970 ¢3 0.691740 -068 B.5411 —0.90 + 0.14 
Flat wow,CDM 0.40370.038 0.5977 9.028 73.3411 —0.93 + 0.15 —0.1735 
External Probes (No SHOES): Flat wCDM 
Planck and Pantheon+ 0.32540 -010 0.675+0-008 66.49+6-39 — 0.9821 0.022 
Planck and galaxy BAO and Pantheon+ 0.31970 -098 0.6817 9:32 66,7818%8 —0.97419 537 
Planck and all BAO and Pantheon+ 0.316 0-605 0.68440 008 66.874 499 —0.978 70 024 
External Probes (No SHOES): Flat wowgCDM 
Planck and Pantheon+ 0.31846.314 0.682+6.014 67.4th4 — 0.8514 8.05 —0.7010-49 
Planck and galaxy BAO and Pantheon+ 0.31875 O08 0.6827 9-008 67.127 9.25 0.878) 5¢5 —0.4543 33 
Planck and all BAO and Pantheon-+ 0.31670 bos 0.684 40.005 67.417 933 —0.8417-066 —0.651935 


Note. Summary of marginalized parameter constraints for Pantheon+ and other external probes. The mean and 68% confidence limit are provided for each 
cosmological parameter. A blank value indicates a parameter is not used in the cosmological fit. 
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Figure 7. Visualizing the impact of a number of the top systematic 
uncertainties in this analysis. The ju residuals are described by Equation (6). 
Each of these systematics is explained in Section 3, and they are combined to 
form the covariance matrix shown in Figure 6. Fragilistic provides nine 
systematic sets of trained SALT2 models, zero-point solutions, and filter 
central wavelengths. Here we show the impact on distance of just the first three. 


and find go=—0.51 + 0.024. In a flat universe gy =“ — 1, 


which gives 2047= 0.326 + 0.016, consistent with the result for 
Qyy reported in this work. Results for Hp from the inclusion of the 
SHOES Cepheid host distances are discussed below. 

The constraints on Qy and Q, for a ACDM model are 
shown in Figure 8. We find (4,;=0.306+0.057 and 
Q, =0.625 + 0.084; a flat universe is within the 68% 
confidence region, and Qy=0 and (Q,=0 are together 
rejected at 4.40 using only the SNe. 

For a flat wCDM model, from the SNelIa alone (not 
including SHOES Cepheid calibration), we find Qy= 
0.3097) :083 and w = —0.90 + 0.14, as shown in the third row of 
Table 3 and in the blue contour of Figure 9. This result is 
consistent within lo of the cosmological constant (w= —1). 
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Figure 8. Confidence contours at the 68% and 95% levels for the Qy, and OQ, 
cosmological parameters for the ACDM from the Pantheon-+ data set, as well 
as from the Planck and combined BAO data sets. The constraints from 
including both the statistical and systematic uncertainties (shaded red) are 
shown as well as when only statistical uncertainties are propagated (unfilled 
dashed). We include two lines for reference: one for a flat universe, where 
Qy + OQ, = 1 and the other that indicates an accelerating universe. 


For a flat wow,CDM model, from the SNe Ia alone (not 
including SHOES Cepheid ano, we find wo= 
—0.93 + 0.15 and w, = —0.1*9, as shown in the fourth row of 
Table 3 and in Figure 10. These results are again consistent 
with a cosmological constant. 

Using distances and a stat-++syst covariance matrix that extends 
to the Cepheid calibrators (Equation (15)) and combining the 
Pantheon+ SNe with the SHOES Cepheid host distance 
calibration, we are able to robustly and simultaneously constrain 
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Figure 9. Sixty-eight percent and 95% confidence contours for flat wCDM for cosmological parameters Q.y, Ho, and w. The contours from the Pantheon+ (red), 
Pantheon+ and SHOES combined data set (teal), Planck Collaboration et al. (2020) TTTEEE-lowE constraints (gray) are shown. The combination of Planck and 
Pantheon-+ (blue) is also shown, which is consistent with a cosmological constant. Planck constraints are bounded by 0.2 < Qy < 0.4 for computational speed. The 


histograms depict marginalized relative probabilities between probes. 


Ho and other cosmological parameters describing the expansion 
history. While we use SHOES Cepheid data and covariance in this 
work, likewise, Pantheon-++ distances and covariance are used 
in R22 in order to fit Ho and qo in flat ACDM. As shown in the 
top Pantheon+ and SHOES section of Table 3, for ACDM, flat 
wCDM, and flat wow,CDM, we find Hp = 73.4 + 1.1, 73.5 + 1.1, 
and 73.34 1.1 kms’ Mpc ', respectively. We note that more 
complex models do not result in decreased Ho constraining power 
from the SNe Ia + Cepheids, while this is not necessarily true for 
other cosmological probes (Section 4.4). 


4.4. Constraints on Cosmological Parameters from Multiple 
Probes 


In this work we combine the Pantheon+ SNe with external 
cosmological probes: CMB from Planck (Collaboration et al. 
2020) TTTEEE-lowE and BAOs from SDSS Main Galaxy 
Sample (Ross et al. 2015), SDSS Baryon Oscillation Spectro- 
scopic Survey (BOSS; Alam et al. 2017), SDSS Extended Baryon 
Oscillation Spectroscopic Survey (eBOSS) Luminous Red Galaxy 
(Bautista et al. 2020), SDSS eBOSS Emission Line Galaxies 
(Bautista et al. 2020), SDSS eBOSS Quasar (Hou et al. 2020), 
SDSS eBOSS Lya (du Mas des Bourboux et al. 2020), all of 
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which have been implemented in CosmoSIS. The aforementioned 
BAO constraints are denoted “all BAO”; we also provide 
constraints from the combination of a spectroscopic redshift 
galaxy-only subset of BAO probes denoted “galaxy BAO.” We 
report constraints in Table 3 for combinations of data sets that are 
deemed compatible and discussed below. 

For a flat wCDM model when combining Pantheon+ and 
Planck, we find w =—0.98210-0%2 and Qy = 0.325700. and 
when further including all BAO, we find w = —0.97819-92 and 
Quy = 0.3161 5-00°” both of which are consistent with the 
cosmological constant at ~3% (Figure 11). As can be seen in 
Figure 9, we do not include SHOES in combinations with 
Planck because these measurements are incompatible (R22). 

For a flat wow,CDM model when combining Pantheon+ and 


Planck, we find wo = —0.8511:995 and w,= —0.70704°, and 
when combining Pantheon+, Planck, and BAO, we find wo = 
—0.84170-966 and wa= = 657338. which is moderately con- 


sistent (20) with a cosmological constant (Figure 12). We note 
that this result is not driven by any single probe. In Figure 10 we 
show constraints for Planck alone and for the combination of 
Planck & Pantheon+. While the broader model freedom of the 
flat wow,CDM allows the Planck-alone Ho to be consistent with 
73kms ! Mpc! owing to degeneracy between Hp and w, (see 
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Table 4 

Sources of Uncertainty 
Description Baseline Systematic (S,,) Oy Cwsys Owsys/Fwstat Awsys 
All Systematics 0.019 0.79 —0.009 
Calibration 
SALT2 Train*ZPT Fragilistic best fit 10 covariance realizations 1/3 each 0.009 0.38 0.000 
SALT2 Method SALT2-B22 JLA SALT2 Surface 1/3 0.008 0.33 0.003 
CSP Tertiary Stars Krisciunas et al. (2017b) Stritzinger et al. (2018) 1 0.003 0.13 —0.003 
HST Calspec 2020 update 5 mmag/7000 A 3 0.003 0.13 —0.006 
Redshifts 
Vpec Map 2M++ 2M++ iLOS & 2MRS 0.7 each 0.002 0.08 0.005 
Redshift Bias No z-shift 10~* z-shift 1 0.011 0.46 0.015 
Astrophysics 
Intrinsic Variations BS21 dust model G10 and C11 0.7 each 0.002 0.08 —0.003 
MW E(B — V) Schlafly & Finkbeiner (2011) 4% Scaling 1.0 0.008 0.33 —0.010 
MW Color Law Fitzpatrick (1999) Cardelli et al. (1989) 1/3 0.006 0.25 —0.006 
Mass Step Split at 10 Split at 10.2 1 0.001 0.04 0.000 
Modeling 
Selection Efficiency Nominal exposure time 20% increase 1 0.004 0.17 0.001 
Populations BS21 parameters Three Sets of Params (P21) 0.6 0.000 0.00 0.003 


Notes. A summary of the systematic uncertainties and the baseline component of the analysis as described in Section 3, the size of the systematic Sy used to determine 
the impact of that systematic, the scaling of the systematic o,, as constrained in this analysis, and the contribution to the total uncertainty in wCDM (can be compared 
to statistical uncertainty of 0.03), and the shift when allowing the uncertainty on the best-fit cosmological parameter. The last column shows the simplistic change in 
best-fit cosmology if a perturbation of size 0, is applied with statistical-only uncertainties. The amount shown is different than seen for the combined shift for the best 
fit and increase of uncertainty given in the previous columns due to the self-calibration, as explained by Brout et al. (2021). 


* ZPT denotes light-curve fitting zero-points. 
Constraints are combined with Planck Collaboration et al. (2020). 


Figure 10), after combining Planck with Pantheon+, the Ho/w, 
degeneracy is broken (Hp = 67.4*}'5 kms 'Mpc_'). Therefore, 
the inclusion of SHOES with Planck & Pantheon-+ results in a 
Bayesian evidence ratio of —9, and we deem this set of probes 
incompatible and do not include them in Figure 10 nor in Table 3. 


4.5. Impact of Systematics on Cosmological Parameter Fits 


To understand the impact of systematic uncertainties, in Table 4 
we group the systematics investigated in this work into four main 
categories: calibration/SALT2, redshifts, astrophysics, and mod- 
eling. The baseline, systematic treatments (S,,), and scaling priors 
(a; as described in detail in Section 3) are summarized for each 
source. The final three columns of Table 4 relate to fits of the 
sample when combined with Planck Collaboration et al. (2020) in 
a flat wCDM model when isolating that systematic. We define 
both the change in the best fit (Aw,,,) and the systematic 
uncertainty contribution to w (037°) as follows: 


Ww 


AWeys = Weys — Wstat (16) 


syS __ 2 — we 
oy = O wrot A wstat > 


(17) 


where Wey, and Oyto are the cosmological constraints when 
utilizing Cotat+sys, and where wera aNd Oy,star are the statistical- 
only constraints when utilizing Cgtat. 

We find that the final systematic uncertainty in w 
(sys = 9.019) is comparable to yet smaller (~80%) than the 
statistical uncertainty, suggesting that the measurement is not 
systematics dominated. The largest contribution to the 
systematic error budget (0.011) is due to the potential for 
redshift-measurement bias. This is followed by the uncertain- 
ties in the Fragilistic calibration offsets and the resulting 
propagation to SALT2 model-training uncertainties and light- 
curve fitting uncertainties (0.009). Additionally important is the 
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conservative uncertainty that was applied owing to the usage of 
the new SALT2 training methodology (0.008) as well as the 
uncertainty in the MW extinction maps (0.008). 

Interestingly, numerous systematic uncertainties are found to 
be negligible (e.g., BS21 parameters, G10 versus C11) in the 
cosmological parameter budget. While certain systematics 
cause redshift-dependent trends as shown in Figure 7, they 
also change the relative scatter of the Hubble residuals. This 
can most easily be seen for the cosmological likelihood values 
(£) for the distances with different intrinsic scatter models 
shown in Table 2. If the baseline analysis is significantly 
preferred (larger £) by the data over one of the analysis 
variants, the impact of that systematic on cosmological 
constraints will be reduced, as is the case for intrinsic scatter. 

As we have built a covariance matrix that includes the 
Cepheid calibrators, we can measure Hy with and without 
systematic uncertainties. For flat ACDM, we find Hy= 
73.6 + 1.1kms 'Mpc™', and when considering only statis- 
tical uncertainties from the SNe alone (excluding Cepheid and 


physical distance calibration uncertainties), on = 0.7 


kms! Mpe™', and oF, = 0.29 km s_'Mpc'. This suggests 
that SN systematic uncertainties are not dominating the 
constraint on Hy and cannot explain the ~7kms ' Mpc! 
difference between Planck and SHOES. 

In Figure 13 we show deviations to the best-fit Hp for each 
individual source of systematic uncertainty relative to the baseline 
analysis and assuming ACDM. For reference, we also show the 
full SN contribution to the Ho error bar (dashed). The deviations 
from the baseline (AH) are small and added in quadrature to 
0.32 kms ' Mpc '. We note that when assessing redshift-specific 
systematics, because model redshifts are not used for the SN- 
Cepheid calibration in Equation (14), they mainly impact the 
Hubble-flow SNe (third rung of the distance ladder). 
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Figure 13. The impact on recovery of Hp, as explained in Section 2, of the 
systematic uncertainties described in Table 4. The units of these measurements 
are kms 'Mpc !. The dashed lines are given at a AHp of 0.7, which is the 
entire contribution of the uncertainty in R22 from SN measurements. 


Finally, to help visualize the impact of systematic uncertain- 
ties, we show in Figure 8 the constraints when including either 
statistical-only uncertainties or the combined statistical and 
systematic uncertainties. Error budgets for different cosmolo- 
gical parameterizations can be generated with the delineated 
files for systematics provided as part of this release. 


4.6. Local Structure in the SN Ia Hubble Diagram 


Large compilations of SN distances have provided impetus 
for searches of local structure, over/underdensities, and proper 
motion (e.g., Mathews et al. 2016; Soltis et al. 2019; Hu et al. 
2020). As an initial study, we create sky maps of the SN 
Hubble diagram residuals (see Figure 14) and examine two 
specific areas on the sky that have been documented in the 
literature and have sufficient SN statistics in the Pantheon+ 
sample for study. 


4.6.1. The CMB Kinematic Dipole 


The motion of the Milky Way and solar system relative to 
the CMB rest frame (v= 369.82 kms_') is corrected for 
following Carr et al. (2021) and Peterson et al. (2022). The 
effect of the CMB dipole motion can be seen in the zyp, sky 
map (middle-right panel of Figure 14), where zp, is the 
heliocentric redshifts. The z¢yp skymap (middle-left panel of 
Figure 14) has the CMB dipole-causing peculiar redshift 
removed, following Equation (7) of Peterson et al. (2022). The 
direction of the CMB dipole, /= 264° and b = 48° (red “o” in 
Figure 14), is shown for reference as well as its antipole 
(red “x’’). 

As discussed in Section 3.1.3, we examine different velocity 
reconstructions due to local structure that include estimates of 
the bulk flow; these are the 2M++ (Carrick et al. 2015) and 
2MRS (Lilow & Nusser 2021) corrections and are shown in the 
top row of Figure 14. These corrections also include the CMB 
dipole correction. Peterson et al. (2022) showed that the 
peculiar-velocity corrections overall reduce the Hubble residual 
scatter by ~ 10%, and this is qualitatively confirmed in our 
maps. The heliocentric map shows a strong dipole as expected; 
the z¢mp Map shows the dipole somewhat removed but with an 
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overcorrection (as expected at low-z because local galaxies 
share some of our motion); and both zyp maps show that the 
peculiar-velocity corrections have removed most of the 
overcorrection. 

However, both reconstructions produce a small signal that 
can be seen in the maps in the direction opposite the motion 
causing the CMB dipole. This signal is found to be local, at 
z< 0.02, and grows with decreasing redshift until z+ 0.01 
(bottom-left panel of Figure 14). A possible reason that there is 
a residual signal in the negative dipole direction in both the 
Zcmp and peculiar-velocity corrected redshifts is that the MW 
motion is coupled with the motion of nearby galaxies in a way 
that is not yet sufficiently modeled. It is also likely that this is 
due to low-number statistics (this is only a lo deviation) and 
the uneven sky coverage (the SNe in this region are mostly 
clustered in Stripe-82). Lastly we note that the positive 
residuals are driven by SNe at z< 0.02, and thus are not 
included in the SHOES (Riess et al. 2022) sample and inference 
of H 0: 


4.6.2. The CMB Cold Spot 


The “CMB cold spot,” a 5° region of —70 zK centered 
at (I~ 209°, b~ —57°), was first detected in data from the 
Wilkinson Microwave Anisotropy Probe (Vielva et al. 2004; 
Cruz et al. 2006), and subsequently in Planck data (Gurzadyan 
et al. 2014). Evidence for an underdensity aligned with the 
CMB cold spot was presented by Rudnick et al. (2007). 
Szapudi et al. (2015) and Kovacs et al. (2021) subsequently 
found the Eridanus supervoid in the direction of the cold spot at 
z 0.15. However, it is not clear if the alignment of Eridanus 
and the CMB cold spot is causal or coincidental. 

We find a signal in the Pantheon+ Hubble diagram when 
examining SNe within a 20° radius of the location of the CMB 
cold spot (blue circle region in the top-left panel of Figure 14). 
The difference in Hubble diagram residuals as a function of 
redshift is shown in the bottom-right panel of Figure 14. There 
are nine SNe in this region of the sky with redshifts on the near 
side (0.12 <z<0.15), and there are 14 SNe on the far side 
(0.15 < z< 0.20) of the proposed void at z= 0.15. There is a 
Hubble residual difference of —0.15 + 0.06 mag between these 
two sets of SNe. For an estimate of the significance, we 
examine 1000 randomly selected 20° apertures across the sky 
with at least eight SNe in each of the near and far redshift 
ranges split on redshifts between 0.08 and 0.20, and find that 
deviations with a similar significance occur only 0.2% of the 
time. We note however, that there are not many independent 
regions that satisfy the selection criteria, and the vast majority 
of the SNe in the cold-spot selection come from the small deep- 
field patch within that region. Taking 100 random samples of 
10° radius from the largest densely sampled region in Pantheon 
+ (Stripe-82 region), we find no other patch has a significance 
that exceeds 1.60, making the Eridanus patch the most 
significant step at that redshift in our data. 


5. Discussion 


This analysis is the latest in a series of papers that attempts to 
both grow the compilation of measured SN Ia light curves and 
improve on the systematic floor. The two most recent 
compilations and analyses are those of JLA and Pantheon, 
which, respectively, included ~ 40% and ~ 60% of the SN light 
curves analyzed here. As seen in Figure | of Scolnic et al. (2022), 
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Figure 14. Healpix (NSIDE = 16) Hubble residual sky maps (the color bar is residual magnitudes) with 20° 2D-Gaussian kernel smoothing, and Hubble residuals for 
two selected apertures. z > 0.01 is applied. Dots show the locations of the SNe in the Pantheon+ sample, with white dots showing the nearby SNe (z < 0.15) and 
black dots showing the distant SNe (z > 0.15). Top left: Hubble diagram corresponding to the baseline analysis utilizing both zcyp dipole corrections and 2M+-+ 
peculiar-velocity corrections. The circled regions designate the 20° regions centered on the negative CMB dipole (red) and CMB cold-spot directions (blue). The small 
circle in the top right (and “x” in bottom left) of each panel represents the direction (and opposite direction) of the motion causing the CMB dipole. Top right: same as 
top left, but instead using 2MRS peculiar-velocity corrections. Middle left) same as top left, but instead not applying any peculiar-velocity corrections. Middle right: 
same as top left, but instead applying neither peculiar-velocity corrections nor the CMB dipole correction. Bottom left: 20° region aligned with the (opposite) CMB 
dipole velocity depicting Hubble diagram residuals as a function of redshift. Bottom right: same as bottom left, but with aperture centered at the CMB cold spot 


(1 = 209°, b = 57°), and over a higher redshift range. 


the majority of the statistical increase for Pantheon+ is in the 
addition of numerous low-redshift samples extending down to 
z=0.001. However, the largest differences in the Hubble 
diagram are not solely the result of statistical increase, but rather 
due to improvements in our methodology. 

We show in Figure 15 the difference in inferred distance- 
modulus values (marginalized over M) for the Pantheon+ sample 
relative to the assumptions used in the JLA analysis, for the three 
most significant improvements presented in this work. First is the 
update in the flux cross-calibration to the Fragilistic solution, 
which impacts both the training of the SALT2 model and the 
zero-points used in light-curve fitting. Second is the impact from 
updating the MW extinction curve used in JLA (Cardelli et al. 
1989) to the Fitzpatrick (1999) relation that is used here. Third is 
the change resulting from improved modeling of the SN Ia 
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intrinsic scatter; while in this work we adopt the BS21 model, we 
include the models developed for JLA (G10 and C11) as 
systematics. Each of these changes has been motivated externally 
by previous works (e.g., Schlafly & Finkbeiner 2011; Brout & 
Scolnic 2021; Brout et al. 2022); however, they nonetheless cause 
shifts in dus/dz of ~0.05, or ~0.04 in w. Finally, because all of 
three of these changes have the same sign of d/dz slope, rather 
than canceling each other, when combined in this work, they 
result in a ~0.1 difference in the constraint on w relative to JLA 
(after combining with CMB). 

As discussed by Scolnic et al. (2019), the constraining power 
of large samples of SNe Ia extends beyond inferences of Ho and 
w/Qy. Large compilations of low-z SNeIa enable precision 
measurements of the local growth-of-structure, typically para- 
meterized by fog (e.g., Huterer et al. 2017; Stahl et al. 2021). 
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Figure 16. Constraints on Qy in flat ACDM when the bounds of the redshift 
cman (C11 Scatter Model range of the sample are changed. In the top panel, the minimum redshift is 
varied. The nominal minimum redshift is 0.01 for Pantheon+ cosmology fits 
without SHOES. In the bottom panel, the maximum redshift is varied. The 
nominal maximum redshift is 2.4 for all fits. 
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compared to the value of (Q,, from the full sample. We show a 
similar test in Figure 16 and find relatively stable values of Qu, 
with no signs of the underdensity seen by Colgain (2019). 
Investigation of the main goal of this work (constraints from 
SNelIa alone for a flat wCDM model) results in stat+syst 


uncertainties of *9-°8 and 0.13 for Qy and w, respectively. This 


represents a factor of two improvement in figure of merit over the 
original Pantheon (stat+-syst uncertainties 0.072 and 0.22 for Qy, 


Hsys 


0.01 0.1 iT and w). This cannot be explained solely by statistical improve- 
Badaliitt ments, but rather is also due to a leap in systematics methodology 

vedsnil over the original Pantheon and JLA. As shown by Brout et al. 

Figure 15. Largest differences in analysis compared to Betoule et al. (2014) (2021), cosmology uncertainty budgets are improved by a factor 
and Scolnic et al. (2018b). Top: updating the extinction curve used in the light- of ~1.5 when not binning or smoothing data and covariance. In 


curve fitting from CCM to F99. Middle: updating the SALT2 model, as 


discussed in Brout et al. (2022). Bottom: changing the baseline assumption for Appendix B we discuss and show a binned error budget for 


the intrinsic scatter to the P21, G10, and C11 models. comparison and find a similar factor of 1.5 improvement from this 
choice alone. In examining the unbinned error budget in Table 4, 

Work is ongoing for this measurement using the Pantheon+ it can be seen that several systematics are no longer impacting 
sample (S. Boruah et al. 2022, in preparation), which will include SN Ia cosmology analyses as strongly as had previously been 
validation with simulations as well as propagation of the thought. One such example is the negligible size of the parent 
covariance matrix, which previously would have had a limited population systematic despite including three additional sources of 
effect on og calculations owing to smoothing/binning over scatter model uncertainty, as was also seen by Popovic et al. 
redshift. (2021a). This, as well as the reduction of a number of other 
While in Section 4 we show a Healpix map of Hubble residuals systematics in comparison to their size in binned analyses (also 
across the sky, there are additional and related tests of anisotropy shown in Appendix Table 6), is due to the power of the large data 
that can be performed with these data. Previous analyses of the sets themselves to self-constrain the size of systematic uncertain- 
first Pantheon sample (e.g., Andrade et al. 2018; Brownsberger ties when the systematic itself is not solely degenerate with the 
et al. 2019; Colin et al. 2019; Soltis et al. 2019) typically search cosmological model parameterization. This is especially important 
for radial or hemispherical residuals across the sky. The addition because it brings this work from potentially being dominated by 
of statistics in the low-redshift sample and improved accounting in systematics to rather being dominated by statistical uncertainties. 
Pantheon-+ would particularly strengthen these types of studies. A Furthermore, as shown by Brout et al. (2021), as data sets grow in 
search for matter over/underdensities was performed by Colgéin size, many systematics will continue to shrink without any 
(2019), which varied the minimum and maximum redshift in the additional effort. Lastly, it is important to note that approaches 
original Pantheon sample and redetermined cosmological con- such as the approximate Bayesian computation method given by 
straints. Colgain (2019) found for Pantheon that 4, could be <0 Jennings et al. (2016) will not be able to make use of this self- 
for a low maximum z of ~0.15, though with only ~ 2c difference constraining benefit unless additional parameters are included to 
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allow the data themselves to scale the input sizes of the systematic 
uncertainties (S,,, in Brout et al. 2021). 

While the SN Ia mass step has received much attention in the 
last decade, we find here that its contribution to the error budget 
is exceedingly small. Unlike previous analyses, the mass-step 
treatment in this work is based on an SN color- and dust- 
dependent model (BS21). We find that this more physical model 
results in smaller scatter in the Hubble diagram (Table 2) and 
better xX? relative to cosmological models, which then results in 
smaller systematic uncertainties. We note that properties of SN Ia 
host galaxies other than stellar mass have been seen to correlate 
with SNIa Hubble diagram residuals. Star formation rate, 
specific star formation rate (sSFR), stellar-population age, and 
metallicity have all been shown to correlate to varying degrees 
with the distance-modulus residuals after standardization 
(Lampeitl et al. 2010; Sullivan et al. 2010; Childress et al. 
2013b; Rigault et al. 2013; Rose et al. 2019). For this reason, 
using sSFR values presented by S22, we also examined the size 
of an sSFR step in the subset of the Pantheon+ sample for which 
we have obtained sSFR measurements (z< 0.2). Without 
applying any bias corrections, we find a significant step in sSFR 
(across the median sSFR) of 0.031 + 0.011. However, after 
applying the nominal set of dust and mass-based bias corrections 
(BS21) used in this analysis, we find a step in sSFR of 
0.008 + 0.011, consistent with zero. This is likely due to galaxy 
properties (i.e., stellar mass) being linked to dust properties, and 
that applying a dust-mass correction is accounting for most, if 
not all, of the correlations with sSFR and is also tracing the dust 
distribution. 

Going forward, statistical constraints on w and 2 from SNe 
will improve significantly owing to upcoming data sets from SN 
programs of the DES (D’ Andrea et al. 2018), Zwicky Transient 
Facility (Dhawan et al. 2022), Young Supernova Experiment 
(Jones et al. 2021), Legacy Survey of Space and Time (The LSST 
Dark Energy Science Collaboration et al. 2018; Sanchez et al. 
2022), Nancy Grace Roman Telescope (Hounsell et al. 2018), etc. 
It is likely that these future data sets will improve the statistical 
precision by a factor of 100 (Scolnic et al. 2018a). 

The size of systematic errors on cosmological parameter 
estimates matched the statistical errors for JLA and the original 
Pantheon. Systematic uncertainties in this work have been 
reduced in comparison to Pantheon, and while their impact is 
still significant, it is no longer the dominant component of the 
total uncertainty. With the coming surveys, systematics will 
also likely improve alongside the increase in statistics, as has 
been the case for previous analyses over the last two decades, 
and as expected from the impact of the systematic “self- 
calibration” described in Brout et al. (2021). 

As shown in the systematics error budget Table 4, the 
dominant sources of systematic uncertainty are now from (1) 
the combination of SALT2 training and calibration of surveys, 
(2) potential redshift-measurement biases, and (3) Milky Way 
dust systematics. Fortunately there are paths forward for each 
of these. For survey flux calibration, dedicated programs are 
needed, and there are currently multiple paths underway to 
improve the fundamental calibration of SN Ia samples and how 
they are tied to various other samples (e.g., Regnault et al. 
2015; Stubbs & Brown 2015; Narayan et al. 2019). There is 
also ongoing work (G. Taylor et al. 2022, in preparation) to 
train the SALT2 model with more photometric systems, which 
has already shown promising improvements to systematic 
uncertainties and the ability to constrain the rest-frame U band. 
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The systematic from the redshift-measurement floor has the 
potential to be reduced using improved cosmology fitting 
methodology, although the extent to which the data itself can 
constrain the size of this floor remains unproven. Alternatively, 
future large surveys can use multiple spectroscopic instruments 
and redshifting codes to mitigate potential sources of redshift- 
measurement bias. The Pantheon+ sample is especially 
sensitive to Milky Way dust systematics because of the 
differences in the samples used for low and high redshift. At 
low redshift, to obtain sufficient statistics in a volume limited 
sample, we have used SNe across the sky and with up to 0.2 in 
MWEBYV, whereas the high-redshift surveys have been carried 
out in low extinction regions of the sky (MWEBV < 0.05). 
Future surveys of larger volumes will be able to mitigate this 
with a plethora of both low and high redshift in low MW 
extinction regions on the sky. 

Throughout this work, there are a number of upstream 
components of this analysis that impact downstream analysis 
steps; i., new calibration (or MWEBV maps/color law) 
motivates new SALT2 training, which motivates new fitting of 
the SN parent populations, which motivates new bias corrections. 
The Pippin framework (Hinton & Brout 2020), used extensively 
in this work, was intentionally developed to automate and 
asynchronize this multistep type of analysis; however, it has yet to 
incorporate aspects such as the SALT2 retraining (Taylor et al. 
2021) or population fitting (Popovic et al. 2021a). Likely, this 
framework will need to expand for future analyses. 

There is an alternate approach to obtaining cosmology 
constraints from SNe that has been gaining traction over the 
last decade. Bayesian hierarchical models have been developed 
that utilize bias-corrected observables (Shariff et al. 2016) and 
that incorporate selection effects directly into the model (Rubin 
et al. 2015) or likelihood (Hinton et al. 2019). However, unlike 
BBC in combination with CosmoSIS, these methods have not 
been validated with large realistic simulations. As noted in 
Appendix C, we release, as part of this analysis, 10 realistic 
simulations of the Pantheon+ data set for such validations. 

While constraints on w_ should easily improve with 
upcoming large SN samples, the road to improving constraints 
on Ho is more challenging. There are a limited number of 
SNe Ia that will explode in the near future within a ~40 Mpc 
radius, a constraint due to HST discovery limits of Cepheids. 
At roughly one SN Ia per year, it will take several decades to 
double the current sample of 42 SNe calibrated by SHOES 
Cepheid hosts. Fortunately, we find that the systematics in the 
measurement of Hy from the SNe are at a scale of 
0.3kms~'Mpce', as shown in Figure 13. This is consistent 
with the general finding of Brownsberger et al. (2021), who 
showed how robust Ho is to systematic uncertainties in 
comparison to the relatively calibration-sensitive constraints of 
Wo or Qy. Lastly, there is ongoing work that combines the 
progress used here by Peterson et al. (2022) and applies it to a 
“two-rung” distance-ladder analysis, in which SNe are 
excluded from the distance ladder (Kenworthy et al. 2022). 


6. Conclusion 


This work is the culmination of a number of supporting 
analyses as part of the Pantheon+ effort. In this work, we 
summarize the various inputs and analyses required to combine 
the supporting works and ultimately measure distances and 
cosmological parameters. For the first time, we are able to 
measure the cosmic expansion history and the local distance 
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ladder Ho simultaneously. We combine our results with 
additional external probes. Importantly, we release a number 
of data and analysis products to facilitate reproducing our work 
by the community. This includes a joint covariance of SNe 
used for measurements of Ho and w. 

For our main results, we find Q,=0.334+0.018 in flat 
ACDM from SNe Ia alone. For a flat woCDM model, we 
measure Wo = —0.90+0.14 from SNela alone and wo= 
—0.978"0 33} when combining SNe with constraints on the 
CMB and all BAO; both are consistent with a cosmological- 
constant model of dark energy. We also present the most 
precise measurements to date on the evolution of dark energy 
in a flat wowgCDM universe, and measure w, =—0.1*$} from 
Pantheon+ alone and w, =—0.65'(38 when combining with 
CMB and BAO data. Finally, while nominal constraints on Ho 
are presented in a companion paper by the SHOES team (R22), 
we perform joint constraints of Ho with expansion history and 
find Hjp=73.5+1.1 in flat wCDM, and we show how 
systematic uncertainties in measurements of the SN component 
of the distance ladder cannot account for the current level of the 
“Hubble tension.” 
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Simulations, light-curve fitting, BBC, and cosmology pipe- 
line are managed by PIPPIN (Hinton & Brout 2020). 
Contours and parameter constraints are generated using the 
CHAINCONSUMER package (Hinton 2016). Plots are generated 
with Matplotlib (Hunter 2007). We use astropy (Price-Whelan 
et al. 2018), SciPy (Virtanen et al. 2020), and NumPy 
(Oliphant 2006). Analysis and visualizations provided in part 
by https: //github.com/bap37 /Midwayplotter. 

D.B. thanks his spouse Isabella and their future daughter for 
their support, as the due date is rapidly approaching! 


Appendix A 
Additional Formalism for Distance and Uncertainty 
Estimation 


As shown in BS21, SN Ia scatter has both a color and host- 
mass dependence (increasing scatter) and a redshift dependence 
that arises from selection effects (decreasing scatter). In this 
work we introduce a new method of accounting for the 
uncertainties using the scatter model predictions. We include 
Oscat(Z, Cc, M,) from simulations as an additive uncertainty 
inside Equation (3) rather than the multiplicative uncertainty 
f(@ c, M,) on the computed ope; that has been used in past 
analyses. The Ogcat(z, c, M,.) term is computed from simulations 
that use the choice of scatter model. The BBC process, after 
correcting distances for selection effects, determines the 
magnitude of Ogca(z, c, M,) in each z, c, M, bin by requiring 
that the observed-simulated distance reduced \7 in each bin is 
unity. If the simulations using a model of intrinsic scatter fully 
describe the observed scatter in the data, the uncertainty 
modeling term in Equation (3), Oscar(z, c, M,), will cause Ogray 
to be 0. 

In the case of the decrease in observed scatter at high redshift 
arising from only intrinsically bright/blue events being 
selected at the limits of the telescope (Kessler et al. 2015), 
we instead apply it as a downscaling of f(z, c, M,) of the 
reported measurement uncertainty and set o.ca(z, c, M,) =0. 
Conversely, for bins of z, c, and M, with xX? greater than unity, 
the necessary Ogcar(z, c, M,) is applied, and f is set to 1. The 
resulting f(z, c, M,) and Ogea(z, c, M,) found from the 
simulations are applied to the Pantheon+ data. 

The method and dimensionality for the application of bias 
corrections are dependent on the adopted scatter model. Table 5 
summarizes the differences between the two main methods 
used in this work, the first of which is applied when assuming 
the BS21/P21 scatter model, and the other when assuming 
the G10 or C11 scatter model. The main difference between 
these groups of scatter models, as discussed in Section 3, is 
whether the intrinsic scatter is driven by diversity in the 
reddening ratios Ry of the light curves, which affects the 
application of bias corrections. For both analysis paths, we 
follow the methodology introduced by Popovic et al. (2021b). 
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Table 5 
Distance Bias (and Uncertainty) Estimation for Scatter Models 


BS21/P21 


G10/Cl11 
Dimensionality 
Mass-step correction 7 a fitted parameter 
Intrinsic Scatter Floor Oroor = Omiy 


Selection Effects 


f(z, ©) 


4D (z, x1, c, M,) 
y corrected for within dpia, (7 and dpost Consistent with zero) 
Oroor = Treat (Zis Cis My) + Cae: applied when f(z, c, M,) > 1 
F(z, c, M,) < 1, applied when ee (zi, ¢;, M,) = 0 


Note. The formalism for four-dimensional and seven-dimensional bias corrections is described by Popovic et al. (2021b) that depend on the intrinsic scatter model 
assumed—either G10/C11 or BS21/P21. The statistical and intrinsic scatter uncertainties from Equation (3) are shown here; the other uncertainty components from 


Equation (3) are independent of the scatter model. 


Appendix B 
Binned Systematic Error Budget 


In Table 6 we show a systematic error budget that is nearly 
identical to what was performed in Table 4, except that the data 
set (AD) and covariance matrix (Ctat+syst) are binned in 20 
redshift bins. This error budget is similar to the methodology 
performed in the most recent SN cosmology analyses where 
binned covariance matrices were used (e.g., Pantheon and 
DES3YR; Brout et al. 2019b) and where smoothed data vectors 
and matrices (which were shown to be equivalent to binned) 
were used (JLA). The total systematic error when binning is a 
factor of 1.5 larger (0.029) than when not binning the data 
set (0.019). 

Systematics that improve the most with unbinned 
matrices are those with smaller owanbinned / owpinned, Binned 
analyses collapse valuable information in the Hubble 
diagram down to a single dimension, redshift. We find that 
as expected, the redshift bias systematic does not improve 
much at all. This is because systematics that only exhibit 
redshift dependence are degenerate with cosmological 
model parameters and cannot be self-constrained by the 
data as easily. Systematics that exhibit dependence in other 
parameters (such as SN color) can be drastically reduced in 
SN Ia cosmological parameter error budgets when not 
performing binned analyses. 
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Table 6 

Comparison of Binned and Unbinned Systematic Error Budgets 
Description AowPinned *owuabinned ayginbinned owginned 
All Systematics 0.029 0.019 0.66 
Calibration 
SALT? Train and ZPT? 0.019 0.009 0.47 
SALT2 Method 0.009 0.008 0.88 
CSP Tertiary Stars 0.005 0.003 0.60 
“HST 0.002 0.003 1.50 
Redshifts 
“Vpec Map N/A 0.002 N/A 
Redshift Bias 0.012 0.011 0.92 
Astrophysics 
Intrinsic Variations 0.009 0.002 0.18 
MW E(B — V) 0.012 0.008 0.67 
MW Color Law 0.007 0.006 0.86 
Mass Step 0.001 0.001 1.00 
Modeling 
Selection Efficiency 0.008 0.004 0.50 
Populations 0.011 0.000 0.00 
Notes. 


* Constraints are combined with Planck prior. 
> ZPT denotes light-curve fitting zero-points. 
“Due to implementation methodology of this systematic, it has not been 
performed in the binned case. 

The increase in the “HST” systematic is likely due to noise, as the values are 
very small for both binned and unbinned. 
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. Ten Catalog Level Simulations of Pantheon+ Light- 
curve Fit Parameters; this work 


The following data products that are provided in part by the 5. os a and Peculiar Velocities; from Carr 
full same of Pantheon-+ supp ieee papers are released P ublicly 6. SN Distance Modulii and Redshifts; this work; Carr et al. 
in machine-readable format’” at pantheonplusshOes.github.io 40 
and as part of SNANA and CosmoSIS (where noted) (er) eee teble § 

: 7. SN Distance Covariance; this work 
1. Light-curve Photometry, Redshifts, and Host-galaxy 8. Cepheid Host Distances; from R22 
Properties; from S22 and Carr et al. (2021) 9. Cepheid Host Distance Covariance; from R22 
2. Trained SALT2-B22 Model; from Brout et al. (2022) 10. SN Ia + Cepheid Host Cosmology Likelihood; this work 
3. SALT2 Fit Parameters; from $22 11. SN Cosmology Chains; this work 
Table 7 
Pantheon+ Hubble Diagram 

CID Survey ZHD O2HD ZCMB ZHEL my" hey ae oe c oe x] Ox, mp Omp 
201 1fe LOSS2 0.00122 0.00084 0.00122 0.00082 9.746 1.516 —0.108 0.040 —0.548 0.134 9.584 0.033 
201 Ife SOUSA 0.00122 0.00084 0.00122 0.00082 9.803 1317 —0.033 0.038 —0.380 0.086 9.784 0.035 
2012cg LOSS2 0.00256 0.00084 0.00256 0.00144 11.470 0.782 0.101 0.018 0.492 0.024 11.816 0.024 
2012cg SOUSA 0.00256 0.00084 0.00256 0.00144 11.492 0.799 0.122 0.039 0.713 0.084 11.880 0.036 
1994DRichmond LOWZ 0.00299 0.00084 0.00299 0.00187 11.523 0.881 —0.112 0.026 —1.618 0.050 11.533 0.032 
1981B LOWZ 0.00317 0.00084 0.0035 0.00236 11.542 0.614 —0.005 0.031 —0.445 0.165 11.664 0.034 
2013aa SOUSA 0.00331 0.00085 0.00478 0.00411 11.207 0.594 —0.104 0.054 0.513 0,152 10.891 0.106 
2013aa CSP 0.00331 0.00085 0.00478 0.00411 11.300 0.580 —0.158 0.036 0.633 0.139 10.844 0.100 
2017cbv CSP 0.00331 0.00085 0.00478 0.00411 11.148 0.578 —0.126 0.032 0.617 0.053 10.773 0.094 
2017cbv CNIa0.02 0.00331 0.00085 0.00478 0.00411 11.258 0.578 —0.096 0.035 0.819 0.066 10.914 0.099 


(This table is available in its entirety in machine-readable form.) 


3° Will be made available after publication. 
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40 
oO 


diag 


mgcorr 1 Table 7 is the error on standardized magnitude from the diagonal 


of the covariance matrix. It is for plotting purposes only and not to be used for 
cosmological fits. 
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